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Abstract:  A sub-cell  tensor  based  technique  for 
modeling  dielectric  interfaces  is  introduced  into  a 
(2,4)  FDTD  method.  For  each  cell  containing  an 
interface,  a tensor  based  method  that  enforces 
continuity  conditions  is  used  to  determine  the  fields 
on  both  sides  the  sloped  interface.  These  fields  are 
then  volumetrically  averaged.  The  approach  is  used 
to  calculate  a corrected  field  value  at  each  grid 
point  of  the  large  fourth-order  stencil.  The 
combined  algorithm  is  computationally 
homogeneous,  unlike  most  previous  algorithms  of 
this  type,  and  thus  lends  itself  to  parallel 
processing.  Additionally,  the  method  may  be  used 
with  other  higher-order  stencils.  The  accuracy  is 
tested  using  the  exact  Mie  series  solution  for 
scattering  from  a dielectric  sphere.  It  is  shown  that 
using  the  (2,4)  tensor  method  results  in  -50-70% 
less  error  than  the  (2,4)  standard  Yee  method  in  the 
vicinity  of  a dielectric  sphere. 

Introduction:  The  finite-difference  time  domain 
(FDTD)  method  is  one  of  the  most  widely 
employed  methods  in  computational 
electromagnetics.  As  it  has  been  pointed  out  in 
many  articles,  the  method  has  problems  when  there 
are  curved  boundaries,  which  are  represented  by 
staircases  on  a Cartesian  grid.  If  continuity 
conditions  are  not  properly  maintained  across  these 
curved  interfaces,  inaccuracies  in  the  field 
components  can  occur.  Nadobny  et.  al.  [1] 
developed  a 3D  tensor  method  for  the  treatment  of 
dielectric  interfaces  to  enforce  continuity  of  the 
appropriate  field  components.  Their  paper  was  a 
major  extension  of  the  work  of  Lee  and  Myung  [2] 
and  demonstrated  much  improved  accuracy  for  the 
standard  (2,2)  algorithm. 

In  this  paper  we  adapt  the  tensor  method  for 
use  with  fourth-order  methods.  Fourth  and  other 
higher  order  methods  permit  modeling  on  coarser 
grids.  This  is  important  because  fourth-order 
methods,  although  very  accurate  in  homogeneous 
regions,  generally  present  accuracy  problems  at 
material  boundaries.  One  remedy  for  this  problem 


is  to  employ  a hybrid  formulation  of  (2,4)  FDTD 
and  sub-grid  (2,2)  FDTD  methods  [3],  where  (2,4) 
stands  for  second-order  accurate  in  time  and  fourth- 
order  accurate  in  space.  In  [3]  a coarse  (2,4)  grid  is 
used  in  the  homogeneous  regions  and  a finer  (2,2) 
sub-grid  near  conducting  walls  and  other  structures. 
Another  method  [4]  uses  a large  (2,4)  region  and  a 
buffer  layer  of  (2,2)  cells  between  the  (2,4)  region 
and  the  interfaces. 

In  [5]  an  efficient  higher-order  alternating- 
direction  implicit  (ADI)  finite-difference  time- 
domain  method  for  unconditionally  stable  analysis 
of  curvilinear  electromagnetic  compatibility  (EMC) 
problems  is  presented.  The  method  is  practically 
dispersionless  and  offers  improved  accuracy  for 
curved  boundaries.  Another  paper  [6]  also  discusses 
the  reduction  of  numerical  dispersion  of  the  finite- 
difference  time-domain  method  based  on  a (2,4) 
computational  stencil.  Rather  than  implementing 
the  conventional  approach,  based  on  Taylor 
analysis  for  the  determination  of  the  finite- 
difference  operators,  two  alternative  procedures  that 
result  in  numerical  schemes  with  diverse  wide-band 
behavior  are  proposed.  The  method  is  shown  to 
outperform  the  standard  (2,4)  method. 

The  method  proposed  here  uses  the  same  (2,4) 
algorithm  and  grid  spacing  for  the  homogenous 
regions  and  across  boundaries  as  opposed  to  mixing 
different  accuracy  (second  and  fourth-order) 
algorithms.  This  is  important  for  parallel 
processing,  i.e.  using  the  Message  Passing  Interface 
(MPI),  where  having  a homogeneous  algorithm  is  a 
great  advantage  so  that  each  processor  executes  the 
same  instructions.  It  also  is  an  advantage 
computationally  if  a fourth-order  accurate  method 
can  be  used  to  model  an  electrically  large  structure 
on  a smaller  coarser  grid,  without  any  special  sub- 
gridding.  We  gauge  the  relative  accuracy  of  the 
standard  Yee  fourth-order  and  combined  sub-cell 
tensor  fourth-order  methods  by  comparing 
computed  results  with  the  exact  Mie  series  solution 
for  plane  waves  scattering  from  a dielectric  sphere. 
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Sub-cell  Tensor  Method:  The  differential  form  of 
Maxwell’s  equations  is  given  by: 


VxE: 

VxH 

where 


(1) 

(2) 


D=6*E  . (3) 

In  the  homogeneous  cells,  where  there  are  no 
interfaces,  Eq.  (3)  can  be  used  to  obtain  E. 
However,  in  those  cells  with  interfaces,  boundary 
conditions  must  be  explicitly  satisfied. 

At  a dielectric  interface  these  continuity 
conditions  must  be  maintained  at  the  interface 
between  media  1 and  2: 


For  Yee  cell  faces  cut  by  an  interface  the 
electric  fluxes  through  the  faces  are  broken  into  two 
parts: 

J \DxdS  = elE\s\+e1Els1x, 

J \DydS  = SlElyS\  + £2E2S2 , (7) 

\\dz<hs = £Xe\s\  + e2e2s2  , 

where  the  S’s  are  areas  and  the  superscripts  stand 
for  side  1 or  2.  This  is  illustrated  in  Fig.  1. 

Combining  Eqs.  (6)  with  (7),  the  following 
tensor  relationship  is  obtained  between  the  average 
electric  flux  density  and  electric  field  in  medium  1, 

D = e-E1  (8) 


(£*iEi  — £*2E2)  ■ n - 0 , (4) 

(Continuity  of  the  normal  components  of  D), 

(E1-E2)xn  = 0,  (5) 

(Continuity  of  the  tangential  components  of  E), 

where  n is  the  unit  normal  vector  to  the  interface. 

Eqs.  (4)  and  (5)  can  be  solved  for  E2  in  terms  of  Ej : 

E2  = AEX  (6) 

where  the  elements  of  the  transformation  matrix  A 
are: 


M 
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where  £ is  a 3 by  3 permittivity  tensor  with 
components: 
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and  where  S = S1  + S2 . 
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Fig.  1.  Center  slice  (face)  of  Yee  cell  centered 
around  Ex(ij,  k) . There  are  two  flux  areas 
separated  by  the  dielectric  interface. 


At  each  point  in  the  stencil  for  updating  the  H 
field  components  the  volume  average  of  the  E 
fields  is  used.  The  fourth-order  update  equation 
for  Hz  obtained  by  discretizing  Eq.  1 is: 

H”+U2(i,j,k)  = H”-U2Q,j,k)  + {-x 

24 

[(-£;  (i,  j + 2 ,k)  + 21E"(i,  j + 1,  k) 

- 21 E” (i,  j, k ) + El (/,  j - 1, k))  - (9) 

(-E;(i  + 2 ,j,k)  + 21E”y{i  + 1 ,j,k) 

- 21  El  (i,  j,  k)  + El  (i  - 1,  j,  £))] 


where  r = 


At 

Axuo 


and  for  example, 


£v  = 


(Fi  El+V2E2X) 


V = Vl+V2, 


(10) 


where  Vl  and  V2  are  the  volumes  on  sides  one  and 
two  of  the  interface.  It  should  be  noted  that  Eq.  (9) 
is  identical  in  form  to  the  standard  (2,4)  update 
equation,  the  only  difference  being  that  each  term  is 
replaced  by  the  volume  averaged  field.  There  are 
analogous  update  equations  for  Hx  and  Hy . 

Standard  fourth-order  update  equations 
for  Dx,Dy  and  Dz  may  be  obtained  by 

discretizing  Eq.  2.  For  example: 
D;+\i,j,k)=D;a,j,k)+^x 

[(-Hnx+V2(i,j,  k + 1)  + 21H”x+ll2(i,j,k) 

- 21Hnx+in(i,  j , k- 1)  + H'l+U2(i,  j,k-  2))  - (11) 
(rHl+m(i  + 1 ,j,k)  + 21Hl+V2(i  + l,j,k) 

- 21Hl+V2(i - 1,  j, k ) + H”z+V2(i - 2,  j, *))] 


where  5 = ^ . 

Ax^o 

The  term  Ex  (/,  j , k)  in  Eq.  (9)  represents  the 
volume  averaged  field  in  the  Yee  cell  centered 
on  Ex  (z,  y , k) . For  any  point  in  the  stencil  with  an 
interface  in  that  cell,  Eq.  (10)  is  used  to  correct  for 


the  interface.  This  concept  is  illustrated  in  Fig.  2.  If 
there  is  no  interface,  then 


This  algorithm  amounts  to  using  a corrected  field 
value  at  each  point  in  the  stencil  to  account  for  any 
interfaces  cutting  through  the  stencil  volume  in  an 
arbitrary  way.  If  the  stencil  volume  has  no 
interfaces  the  algorithm  reduces  to  the  standard 
(2,4).  The  entire  algorithm  may  be  briefly 
summarized  as  follows  for  one  update: 

[1]  Perform  standard  (2,4)  update  of  D 
using  Hx  , Hy  and  Hz . 

[2]  Test  all  8 £ cells  for  interfaces  within  the 
fourth-order  stencil  for  updating  H . 

If  E cell  has  an  interface  then  use  the  sub-cell 
tensor  method: 

(a)  Compute  electric  field  from  average 
electric  flux  density,  E{  ~ £-1  D . 

(b)  Obtain  electric  field  on  other  side  of 

interface,  E2  = . 

(c)  Volume  average  electric  field, 

- (v.El  + V2E2X) 

x V 

else  if  the  E cell  has  no  interface  then, 

(a)  Compute  electric  field  from  Ex  = 

e 

[3]  Perform  standard  (2,4)  update  for  H using 
Ex , Ey  and  Ez . 


Fig.  2.  The  stencil  for  the  (2,4)  FDTD  method 
showing  the  8 E field  cells  (squares)  and  H 
field  cell  (circle).  Two  of  the  E field  cells 
are  cut  by  an  interface  (dotted  line)  at  an 
angle  and  require  the  sub-cell  corrections. 
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Computational  Cases:  Six  computational  cases  are 
performed  to  test  the  accuracy  of  the  standard  (2,4) 
and  tensor  (2,4)  algorithms.  Scattering  problems  are 
done  with  a plane  wave  scattering  from  a dielectric 
sphere  with  a large  dielectric  change  to  emphasize 
errors  near  the  interface.  The  incident  plane  wave  is  ^ 
polarized  in  the  z direction  and  travels  in  the  +y  J 
direction,  as  illustrated  in  Fig.  3.  The  total-  g 
field/scattered-field  (TF/SF)  formulation  is  used  to  ^ 
introduce  a plane  wave  into  the  volume.  Uniaxial 
perfectly  matched  layers  (UPML),  10  cells  wide, 
are  used  for  the  absorbing  boundaries.  Case  I uses  a 
sphere  with  a relative  dielectric  constant  of  4 and  a 
uniform  grid  spacing  of  10  points  per  wavelength 
(ppw)  in  the  sphere.  The  parameters  for  various 
cases  are  summarized  in  Table  I. 


Table  I.  Parameters  for  Cases. 


Freq. 

(GHz) 

Dielectric 

Constant 

ppw 

Sphere 

Radius 

Case  I 

5.0 

4 

10 

7.5 

Case  II 

5.0 

4 

20 

15.0 

Case  III 

5.0 

8 

10 

7.5 

Case  IV 

5.0 

8 

20 

15.0 

Case  V 

5.0 

12 

10 

7.5 

Case  VI 

5.0 

12 

20 

15.0 

Z 


Fig.  3.  Diagram  of  incident  wave,  dielectric  sphere, 
coordinate  system  and  y-cut  at  z value. 

Fig.  4 shows  a typical  computed  cut  for  Case 
III,  parallel  to  the  y-axis  and  through  grid  point 
(x=0,z=l),  near  the  sphere  center.  Shown  is 
Ey(0,y,l)  computed  using  the  (2,4)  tensor  method 
and  the  (2,4)  standard  Yee  method  against  the  exact 
Mie  series  solution.  Fig.  4 shows  that  the  (2,4) 
tensor  method  agrees  much  better  with  the  exact 
solution  than  the  (2,4)  standard  method,  along  the 
entire  cut. 


Ey 


Y-Distance 

Fig.  4.  Case  III.  Comparison  of  (2,4)  tensor  and 
(2,4)  standard  methods  with  exact  Mie 
series.  The  sphere  lies  between  grid  points 
24  and  40, 

Fig.  5 shows  a comparable  cut  for  Ey(0,y,5)  for 
Case  V.  The  (2,4)  tensor  method  is  closer  to  the 
exact  solution  inside  the  sphere,  at  the  sphere 
boundaries,  and  outside  the  sphere.  The  standard 
Yee  method  also  exhibits  pronounced  overshoots  at 
the  interfaces. 


Ey 


Fig.  5.  Case  V.  Comparison  of  (2,4)  tensor  and 
(2,4)  standard  methods  with  exact  Mie 
series.  The  sphere  lies  between  grid  points 
56  and  72. 


Error  Evaluation:  In  order  to  assess  the  relative 
errors  of  the  tensor  and  standard  methods  a 
numerical  comparison  is  made  between  the 
computed  values  and  the  exact  Mie  series  solution. 
The  solution  is  computed  in  spherical  coordinates 
and  transformed  into  Cartesian  coordinates  along 
cuts  through  the  sphere  (shown  in  Fig.  3),  to 
correspond  in  space  to  the  FDTD  spatial  cuts. 
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The  following  error  measure  function  is  used: 


error  = 


^ | & exact 


computed 


IK 


(13) 


This  error  function  is  computed  along  a cut  through 
the  sphere  and  extending  1 radius  beyond  the 
sphere  boundary  on  both  sides  so  that  the  cuts  are  4 
radii  in  length.  The  exact  value  is  taken  to  be  the 
average  of  the  analytical  value  at  center  of  the  Yee 
cell,  obtained  from  the  Mie  series  solution.  Using 
the  spatial  average  of  the  exact  solution  is 
necessary  near  the  jump  discontinuities  to  properly 
compare  to  the  computed  values  which  are  really 
the  average  values  at  the  center  of  the  Yee  cells. 
Table  II  shows  the  errors  computed  for  the  six  cases 
along  y-cuts  at  4 different  z values  for  increasing 
dielectric  constants.  The  ratio  r shown  is  the  tensor 
average  error  divided  by  the  standard  average. 
Table  II  shows  that  r is  about  0.3  at  20  ppw  and  .6 
at  10  ppw  for  the  six  cases.  Cases  II  and  IV  both 
show  a decrease  of  about  50%  in  the  average  error 
by  going  from  10  ppw  to  20  ppw  for  the  tensor 
method.  The  standard  method  shows  worse 
convergence.  Case  VI  shows  only  about  a 20% 
decrease  in  the  average  error  by  going  to  20  ppw. 


Table  II.  Errors  for  Cases  I - VI. 


Case  I 

ez(l) 

ez(4.5) 

ey(l) 

ey(4) 

av. 

ten(2,4) 

1.79 

2.30 

2.31 

2.11 

2.13 

sta(2,4) 

1.94 

3.53 

4.85 

4.06 

3.60 

r=.59 

Case  II 

ez(l) 

ez(8.5) 

ey(l) 

ey(8) 

av. 

ten(2,4) 

0.61 

0.88 

1.46 

0.98 

0.98 

sta(2,4) 

1.80 

1.78 

4.83 

4.41 

3.21 

i=.31 

Case  HI 

ez(l) 

ez(4.5) 

ey(l) 

ey(4) 

av. 

ten(2,4) 

3.98 

5.62 

2.54 

3.13 

3.82 

sta(2,4) 

6.67 

11.30 

6.99 

5.68 

7.66 

i=.50 

Case  TV 

ez(l) 

ez(8.5) 

ey(i) 

ey(8) 

av. 

ten(2,4) 

2.15 

2.16  1 

1.94 

1.48 

1.93 

sta(2,4) 

4.64 

4.95 

6.26 

7.00 

5.71 

i=.34 

Case  V 

ez(l) 

ez(4.5) 

ey(l) 

ey(4) 

av. 

ten(2,4) 

3.45 

3.25 

2.77 

3.26 

3.18 

sta(2,4) 

3.86 

5.33 

8.31 

5.32 

5.71 

r=.56 

Case  VI 

_ez(l) 

ez(8.5) 

_ ey(!) 

ey(8) 

av. 

ten(2,4) 

2.66 

2.90 

2.67 

1.67 

2.48 

sta(2,4) 

5.87 

3.55 

8.48 

8.88 

6.70 

i=37 

Efficiency:  The  (2,4)  tensor  is  compared  with  the 

(2,4)  standard  Yee  for  total  CPU  time  and 
additional  memory  requirements.  The  computations 
were  all  performed  on  an  IBM  p690  parallel 
computer  using  16  processors.  This  case  uses  a grid 
size  of  168x168x168  and  2000  time  steps.  The 

(2,4)  tensor  method  uses  373  s of  CPU  time 
compared  to  251s  for  the  (2,4)  standard  Yee,  or  a 
ratio  of  1.49.  The  tensor  method  also  requires  some 
additional  memory  primarily  to  store  the  9 £~l 
tensor  components  and  the  9 transformation  matrix 
components  of  the  A matrix  for  each  interface  cell. 
For  this  case  there  are  786  interface  cells  for  each 
of  the  staggered  field  positions  Ex , Ey  and  Ez . 
Also  the  volume  fractions  must  be  stored  for  each 
interface  cell.  The  total  additional  memory 
overhead  compared  with  the  standard  (2,4)  method 
amounts  to  only  about  0.2  Mbytes  for  this  case. 

Conclusions:  A tensor  method  to  handle  dielectric 
interfaces  has  been  combined  in  a straightforward 
way  with  a standard  (2,4)  FDTD  algorithm  and 
results  in  a computationally  homogeneous 
algorithm  suitable  for  parallel  computing.  The 
numerical  cases,  using  scattering  from  a dielectric 
sphere,  demonstrate  that  the  combined  (2,4)  tensor 
method  significantly  improves  the  accuracy  of  the 

(2,4)  standard  Yee  method  near  interfaces.  The 
tensor  method  may  be  combined  with  any  higher- 
order  FDTD  algorithm,  involving  a large  stencil. 
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